LA Street Slopes

<<<<<<< HEAD

Project Background

=======

Project Background

>>>>>>> ab384e6cd3f017dce9e7d273fac41bf48f398611

For the LADOT Play Streets Program, I was asked if we had a layer that had the following information:

1. Street Classification
2. Street Width
3. Street Slope

The street classification system is part of the street centerline file, maintained by the Los Angeles Department of Public Works Bureau of Engineering and available at the City’s GeoHub here. Street width is contained in a slightly different version of the same centerline file, maintained by the Bureau of Street Services, and available as part of the ‘Street Pavement Condition’ layer on the City’s GeoHub here.

However, there is not readily available layer with street slope already calculated, so I decided to create one myself, using this blog post as a guide. The author, Stephen Von Worley, combined the National Elevation Dataset’s 1/3-second data with the OpenStreetMap grid to produce his result. For this project, I combined the same elevation data with the two layers I mentioned above. Another option in R that I could have pursued included utiling the R package elevatr, which includes querys for elevation data for a bunch of different APIs. Most of the APIs don’t provide access to Digital Elevation Models, which is the most precise elevation data, and definitely necessary for this project. Mapzen, though it does provide API access to raster layers and DEM-based querying API, it has rate limits, so in this case it was just easier to download the DEM raste data and query from there.

<<<<<<< HEAD

Before starting this project, I went ahead and downloaded the two street centerline shapefiles as well as the “USGS NED 1/3 arc-second 2013 1 x 1 degree ArcGrid” elevation data. For the City of Los Angeles, I needed to download two data elevation tiles: n35w119 (Northern LA), and n34w119 (Southern LA). These NED data are at resolutions of 1/3 arc-second (approx. 10 meters) and in some areas at 1/9 arc-second (approx. 3 feet). This is precise enough for the task at hand.

Project Plan

In order to calculate the slope for each street segment, my plan was to calculate the elevation at each end point of the street segment. I could generate an endpoint file from the centerline network, but there is already one available on the GeoHub here. This intersections layer has IDs that can be easily joined to the street centerline file for the slope calculation.

I originally wanted to use the BSS PCI layer features and join the relevant attributes of the BOE centerline layer to the BSS features. After doing the first exercise, however, I realized that the two layers weren’t quite as compatable as I thought. My linking ID, ‘SECT_ID’ is actually a 1:many join with the BOE centerline – there are several cases where there is the same ‘SECT_ID’ for multiple line features. Unfortunately, to my knowledge, there is not an ID that provides a 1:1 join between the BOE and BSS centerline layers. After seeing some confusing results, I decided I needed to pick one to primarily rely on. If I did the BSS, I would need to calculate the endpoints of each segment, and then derive the elevation from there. If I went with BOE, I would need to calculate the length of each segment. I decided to go with BOE as the primary centerline file, and then just pull the BSS file for the street width later.

=======

Before starting this project, I went ahead and downloaded the two street centerline shapefiles as well as the “USGS NED 1/3 arc-second 2013 1 x 1 degree ArcGrid” elevation data. For the City of Los Angeles, I needed to download two data elevation tiles: n35w119 (Northern LA), and n34w119 (Southern LA). Also, in order to calculate the slope for each street segment, my plan was to calculate the elevation at each end point of the street segment. I could generate an endpoint file from the centerline network, but there is already one available on the GeoHub here. This intersections layer has IDs that can be easily joined to the street centerline file for the slope calculation.

>>>>>>> ab384e6cd3f017dce9e7d273fac41bf48f398611

Project Tools

• sf support for simple features; here is a good overview of the package

• raster For the analysis of the elevation data

• tidyverse Collection of R packages for datascience, including tidyr and dplyr

• leaflet for viewing our output on a web map

<<<<<<< HEAD

## Linking to GEOS 3.6.1, GDAL 2.2.0, proj.4 4.9.3

## Warning: package 'dplyr' was built under R version 3.4.2

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## Warning: package 'geosphere' was built under R version 3.4.2

=======

# load libraries
library(sf)

## Linking to GEOS 3.6.1, GDAL 2.2.0, proj.4 4.9.3

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 3.4.2

## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr

## Conflicts with tidy packages ----------------------------------------------

## filter(): dplyr, stats
## lag():    dplyr, stats

library(leaflet)

# import street data
bss <- st_read("Data/Street_Pavement_Condition.shp",quiet=TRUE)
boe <- st_read("C:/Users/dotcid034/Documents/NonGitHubData/boe/Streets_Centerline.shp", quiet=TRUE)
int <- st_read("Data/Intersections.shp", quiet=TRUE)

>>>>>>> ab384e6cd3f017dce9e7d273fac41bf48f398611

Explore the Data

BOE Centerline

<<<<<<< HEAD

## Observations: 84,186
=======
glimpse(boe)
## Observations: 84,184
>>>>>>> ab384e6cd3f017dce9e7d273fac41bf48f398611
## Variables: 39
## $ OBJECTID   <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, ...
## $ ASSETID    <dbl> 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 1...
## $ INT_ID_FRO <dbl> 4176, 4299, 4373, 4465, 4078, 4149, 4479, 4883, 489...
## $ INT_ID_TO  <dbl> 4259, 1843, 4390, 4488, 4117, 4159, 4624, 4890, 489...
## $ STNUM      <dbl> 1736, 959, 7854, 3463, 2412, 1004, 8609, 6240, 7935...
## $ MAPSHEET   <fctr> 099B209, 088-5A199, 120B177, 117B181, 084B209, 099...
## $ ID         <dbl> 7392, 7417, 7460, 7585, 7015, 7049, 7762, 8129, 813...
## $ ADRF       <dbl> 8301, 701, 2901, 3601, 11501, 1001, 5965, 4601, 440...
## $ ADRT       <dbl> 8399, 799, 2929, 3799, 11599, 1199, 5999, 4699, 479...
## $ ZIP_R      <dbl> 90001, 90044, 90016, 90016, 90059, 90044, 90016, 90...
## $ ADLF       <dbl> 8300, 700, 2900, 3600, 11500, 1000, 5964, 4600, 440...
## $ ADLT       <dbl> 8398, 798, 2928, 3798, 11598, 1198, 5998, 4698, 479...
## $ ZIP_L      <dbl> 90001, 90044, 90016, 90016, 90059, 90044, 90016, 90...
## $ TDIR       <fctr> S, W, S, S, S, W, W, S, S, S, W, S, W, W, S, S, S,...
## $ STNAME     <fctr> WADSWORTH, 110TH, COCHRAN, POTOMAC, SUCCESS, 82ND,...
## $ STSFX      <fctr> AVE, ST, AVE, AVE, AVE, ST, ST, PL, ST, AVE, ST, A...
## $ SFXDIR     <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ STNAME_A   <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ STSFX_A    <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ STATUS     <fctr> O, O, O, O, O, O, O, O, O, O, O, O, O, O, O, O, O,...
## $ TEMP_      <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ SECT_ID    <fctr> 5776100, 6773300, 1229900, 4363900, 5149200, 66710...
## $ DEDRQ      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, ...
## $ REMARKS    <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ SV_STATUS  <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ ST_SUBTYPE <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ CRTN_DT    <fctr> 1995-04-06T00:00:00.000Z, 1995-03-24T00:00:00.000Z...
## $ LST_MODF_D <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ OLD_STREET <fctr> Local Street, Local Street, Local Street, Local St...
## $ PLANNING_S <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ DEDRQ_DESC <fctr> No, No, No, No, No, No, No, No, No, No, No, Yes, N...
## $ TOOLTIP    <fctr> WADSWORTH AVE\nStreet Designation: Local Street - ...
## $ NLA_URL    <fctr> navigatela/reports/centerline_mb.cfm?pk=168&t=5776...
## $ Planning_A <dbl> 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 1...
## $ TYPE       <dbl> 70, 70, 70, 70, 70, 70, 70, 60, 70, 70, 60, 40, 70,...
## $ MODIFIED   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ...
## $ Street_Des <fctr> Local Street - Standard, Local Street - Standard, ...
## $ Street_D_1 <fctr> Local Street - Standard, Local Street - Standard, ...
<<<<<<< HEAD
## $ geometry   <simple_feature> MULTILINESTRING ((-118.2585..., MULTILIN...
BSS Centerline / Street Pavement Condition
=======
## $ geometry   <simple_feature> MULTILINESTRING((-118.25858..., MULTILIN...

BSS Centerline / Street Pavement Condition

glimpse(bss)

>>>>>>> ab384e6cd3f017dce9e7d273fac41bf48f398611

## Observations: 72,113
## Variables: 18
## $ FID        <dbl> 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 100...
## $ SECT_ID    <fctr> 5542400, 5541300, 5541100, 3792300, 4688700, 32909...
## $ STREET_DES <fctr> Secondary Highway, Secondary Highway, Secondary Hi...
## $ PRIME      <fctr> VANOWEN ST, VANOWEN ST, VANOWEN ST, MORSE AVE, ST ...
## $ NAME       <fctr> VANOWEN ST, VANOWEN ST, VANOWEN ST, MORSE AV, SAIN...
## $ OTHER      <fctr> NA, NA, NA, NA, NA, NA, LITTLE, NA, NA, NA, NA, NA...
## $ FROM_      <fctr> WOODMAN AV, GOODLAND AV, BABCOCK AV, VANOWEN ST, V...
## $ TO_        <fctr> MAMMOTH AV, COLDWATER CANYON AV, BELLAIRE AV, WELB...
## $ SURFACE    <fctr> AS, AS, AS, AS, AS, AS, AS, AS, PC, AS, AS, AP, AS...
## $ SURF1      <fctr> Asphalt, Asphalt, Asphalt, Asphalt, Asphalt, Aspha...
## $ STREET_TYP <fctr> SE, SE, SE, LO, LO, LO, SE, LO, LO, LO, LO, SE, LO...
## $ ST_TYPE    <fctr> Select, Select, Select, Local, Local, Local, Selec...
## $ NEWW_CONDI <dbl> 100, 70, 69, 91, 100, 89, 73, 90, 86, 67, 100, 90, ...
## $ STATUS_2   <fctr> Good, Fair, Fair, Good, Good, Good, Good, Good, Go...
## $ ZLENGTH    <dbl> 325, 570, 650, 310, 660, 1955, 1150, 270, 900, 380,...
## $ ZWIDTH     <dbl> 65, 66, 66, 36, 36, 36, 60, 40, 27, 36, 46, 62, 40,...
## $ Shape__Len <dbl> 0.0010913980, 0.0018651708, 0.0021641328, 0.0008528...
<<<<<<< HEAD
## $ geometry   <simple_feature> LINESTRING (-118.4312657132..., LINESTRI...

Data Prep

For the centerline files, both have quite a bit of unneeded information that we can trim out when we do the merge. Since the BSS version of the centerline file contains the street width and the BOE version contains the street classification, we need to merge them (and trim out any other unneeded columns) into one street centerline file. The BSS centerline file is slightly more restrictive – only containing streets that are maintained by the City of Los Angeles – so it will have fewer streets compared to the BOE one. I’ll merge the BOE file with that one.

Note: Once I started bringing in the raster data, I began getting an error while runing the select() function. I found out from this blog post that the error is a result of several packages having a select() function (in this case the function from the raster package is superceding the dplyr one). The solution is simple – make sure to explicitly call out which package you are using with following syntax: package::function()

For the intersection layer, the I’m just going to trim away unneeded columns. First, a glimpse:

======= ## $ geometry <simple_feature> LINESTRING(-118.43126571321..., LINESTRI...

Clean / Prepare Intersection & Centerline Files

For the centerline files, both have quite a bit of unneeded information that we can trim out when we do the merge. Since the BSS version of the centerline file contains the street width and the BOE version contains the street classification, we need to merge them (and trim out any other unneeded columns) into one street centerline file. The BSS centerline file is slightly more restrictive – only containing streets that are maintained by the City of Los Angeles – so it will have fewer streets compared to the BOE one. I’ll merge the BOE file with that one.

Tip: Once I started bringing in the raster data, I began getting an error while runing the select() function. I found out from this blog post that the error is a result of several packages having a select() function (in this case the function from the raster package is superceding the dplyr one). To solve the mixup, just make sure to explicitly designate which package the select() function is from using the syntax package::function() when it is not clear.

# trim bss centerline
bss <- bss %>%
  dplyr::select(SECT_ID, PRIME, FROM_, TO_, SURF1, ZLENGTH, ZWIDTH, geometry) %>%
  rename(sect.id = 'SECT_ID',
         street = 'PRIME',
         from = 'FROM_',
         to = 'TO_',
         surface = 'SURF1',
         length.ft = 'ZLENGTH',
         width.ft = 'ZWIDTH',
         geom = 'geometry') %>%
  mutate(sect.id = as.double(sect.id))

# trim boe centerline
boe <- boe %>%
  dplyr::select(ASSETID, INT_ID_FRO, INT_ID_TO, SECT_ID, OLD_STREET, Street_Des) %>%
  rename(asset.id = 'ASSETID',
         fromint.id = 'INT_ID_FRO',
         toint.id = 'INT_ID_TO',
         sect.id = 'SECT_ID',
         old.desig = 'OLD_STREET',
         new.desig = 'Street_Des') %>%
  mutate(sect.id = as.double(sect.id)) %>%
  st_set_geometry(NULL)

# left join bss <- boe
streets <- bss %>%
  inner_join(boe, by='sect.id')

Let’s take a peek at the merged dataframe:

glimpse(streets)

## Observations: 77,801
## Variables: 13
## $ sect.id    <dbl> 54605, 54594, 54592, 37482, 46329, 32500, 47519, 50...
## $ street     <fctr> VANOWEN ST, VANOWEN ST, VANOWEN ST, MORSE AVE, ST ...
## $ from       <fctr> WOODMAN AV, GOODLAND AV, BABCOCK AV, VANOWEN ST, V...
## $ to         <fctr> MAMMOTH AV, COLDWATER CANYON AV, BELLAIRE AV, WELB...
## $ surface    <fctr> Asphalt, Asphalt, Asphalt, Asphalt, Asphalt, Aspha...
## $ length.ft  <dbl> 325, 570, 650, 310, 660, 1955, 1150, 270, 900, 380,...
## $ width.ft   <dbl> 65, 66, 66, 36, 36, 36, 60, 40, 27, 36, 46, 62, 40,...
## $ asset.id   <dbl> 8537, 13523, 6043, 33706, 49075, 51045, 67151, 5400...
## $ fromint.id <dbl> 24957, 26361, 26333, 37631, 20540, 39573, 13335, 56...
## $ toint.id   <dbl> 24940, 26352, 26329, 37630, 56461, 39396, 13024, 52...
## $ old.desig  <fctr> Secondary Highway, Secondary Highway, Secondary Hi...
## $ new.desig  <fctr> Avenue II, Avenue II, Avenue II, Local Street - St...
## $ geom       <simple_feature> LINESTRING(-118.43126571321..., LINESTRI...

For the intersection layer, the I’m just going to trim away unneeded columns. First, a glimpse:

glimpse(int)

>>>>>>> ab384e6cd3f017dce9e7d273fac41bf48f398611

## Observations: 60,627
## Variables: 17
## $ OBJECTID   <dbl> 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 100...
## $ ASSETID    <dbl> 96967, 96968, 96969, 96970, 96971, 96972, 96973, 96...
## $ CL_NODE_ID <dbl> 13589, 13634, 13660, 15007, 15024, 15030, 15031, 15...
## $ X          <dbl> 6392995, 6406002, 6398905, 6480819, 6479671, 647626...
## $ Y          <dbl> 1839536, 1841805, 1841108, 1861803, 1861156, 186087...
## $ LAT        <dbl> 34.04585, 34.05228, 34.05026, 34.10803, 34.10624, 3...
## $ LON        <dbl> -118.5568, -118.5139, -118.5373, -118.2671, -118.27...
## $ TYPE       <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ CRTN_DT    <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ LST_MODF_D <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ USER_ID    <fctr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ FROM_ST    <fctr> PASEO MIRAMAR, VILLA WOODS PL, ALCIMA AVE, WEST SI...
## $ TO_ST      <fctr> ESPERA AVE, D/E, LAS LOMAS AVE, ROWENA AVE, PANORA...
## $ TOOLTIP    <fctr> PASEO MIRAMAR at ESPERA AVE, VILLA WOODS PL at D/E...
## $ ZIP_CODE   <dbl> 90272, 90272, 90272, 90039, 90039, 90027, 90027, 90...
## $ NLA_URL    <fctr> navigatela/reports/intersection_report.cfm?pk=1358...
<<<<<<< HEAD
## $ geometry   <simple_feature> POINT (-118.556805744008 34..., POINT (-...

Now I’ll trim away the unneeded columns:

int <- int %>%
  dplyr::select(ASSETID, CL_NODE_ID, FROM_ST, TO_ST, X, Y, LAT, LON) %>%
  rename(asset.id = 'ASSETID',
         clnode.id = 'CL_NODE_ID',
         from = 'FROM_ST',
         to = 'TO_ST',
         lat = 'LAT',
         lon = 'LON')

Load Raster Layers

In order to calculate the slope of each street segment, I first need to calculate the elevation at each of these street segment endpoints. In order to get the elevation at each of the endpoints (intersections), I am going to extract this elevation data from the raster elevation image that I pulled from the National Elevation Dataset.

======= ## $ geometry <simple_feature> POINT(-118.556805744008 34...., POINT(-1...

Now I’ll trim away the unneeded columns:

int <- int %>%
  dplyr::select(ASSETID, CL_NODE_ID, X, Y, LAT, LON, geometry) %>%
  rename(asset.id = 'ASSETID',
         clnode.id = 'CL_NODE_ID',
         x = 'X',
         y = 'Y',
         lat = 'LAT',
         lon = 'LON',
         geom = 'geometry')

Extract Elevation Data from Raster Layer

In order to calculate the slope of each street segment, I first need to calculate the elevation at each of these street segment endpoints. In order to get the elevation at each of the endpoints (intersections), I am going to extract this elevation data from the raster elevation image that I pulled from the National Elevation Dataset.

As an aside: Another (perhaps easier) method that I could have used to extract elevation data was through the elevatr package, which

>>>>>>> ab384e6cd3f017dce9e7d273fac41bf48f398611

library(raster)

## Warning: package 'raster' was built under R version 3.4.2

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

<<<<<<< HEAD

# import elevation data & merge into one raster layer
la_n <- raster("C:/Users/Tim/Documents/data/street-slope-la/n35w119/grdn35w119_13")
la_s <- raster("C:/Users/Tim/Documents/data/street-slope-la/n34w119/grdn34w119_13")
la_elev <- merge(la_n, la_s)

# preview the merged raster
plot(la_elev)

Extract Elevation Data

Note: the sf package was released on CRAN on January 2017, while the latest release of the raster package was back in June 2016. Since there is no interoperability between these two packages, I add the extra step of first converting the sf object into an sp object before working with the raster image. According to the documentation, the elevation values are in meters.

# cast as sp object, extract elevation, convert back to sf object
int_sp <- as(int, "Spatial")
int_elev <- extract(la_elev,
                    int_sp,
                    method='simple')

## Warning in .local(x, y, ...): Transforming SpatialPoints to the CRS of the
## Raster

int_elev <- cbind(int, int_elev)
int_elev <- int_elev %>%
  dplyr::select(clnode.id,
                lat,
                lon,
                int_elev) %>%
  st_set_geometry(NULL)

glimpse(int_elev)

## Observations: 60,627
## Variables: 4
## $ clnode.id <dbl> 13589, 13634, 13660, 15007, 15024, 15030, 15031, 150...
## $ lat       <dbl> 34.04585, 34.05228, 34.05026, 34.10803, 34.10624, 34...
## $ lon       <dbl> -118.5568, -118.5139, -118.5373, -118.2671, -118.270...
## $ int_elev  <dbl> 113.78664, 118.57660, 111.74510, 137.62737, 138.3995...

Street Slope

Now that I have the elevation at each of the intersection points, I can calculate the slope of the each street segment by joining the intersection to the (updated) street centerline file. I need to do the join twice – once for each side of the street. Before that though, I am going to calcualte the length of the BOE segment, which also requires projecting it. During this process, I am also going to filter out any of the private / closed / unknown roads.

Steepest streets in LA?

Let’s take a look at the streets in LA with the highest grade. For this exercise, I wanted to look at segments that were at least 150m long.

Winners listed below:

• Mateo: 37.78% This is actually not a true street! It has been permitted for development (which is why the BOE centerline contains the outline for the street), but it doesn’t count for our purposes.

1. Escalon: 33.93% The winner is out in Encino, just above Mulholland Dr.
2. Quito: 30.83% This street exists, and is public, but too narrow for Google car to have gone up.
3. Oletha: 30.51% This is just one street down from Quito in the same Beverly Glen neighborhood.
4. Avenida De Santa Ynez: 29.57% Also on the westside, out near Malibu.

The Webmap

# Set Color Palette
pal <- colorNumeric(
  palette = "YlGnBu",
  domain = streets$slope
)

# Create the leaflet map
map <- leaflet() %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  addPolylines(
    color = ~pal(slope),
    weight = 1,
    opacity = 1,
    data = streets) %>%
  addLegend(pal = pal, 
          values = streets$slope, 
          position = "bottomright", 
          title = "Slope",
          labFormat = labelFormat(suffix = "%"))

map

=======

## The following object is masked from 'package:tidyr':
## 
##     extract

# import elevation data
la_n <- raster("C:/Users/dotcid034/Documents/NonGitHubData/n35w119/grdn35w119_13")
la_s <- raster("C:/Users/dotcid034/Documents/NonGitHubData/n34w119/grdn34w119_13")

-look at documentation for interpreting the extraction

Fun Aside: What are the steepest streets in LA?

-include a google streetview link at the locations

Display the Map

# map <- leaflet() %>%
#   addProviderTiles(providers$Stamen.TonerLite) %>%
#   addPolylines(
#     color = '#0E016F',
#     weight = 1,
#     opacity = 1,
#     data = streets
#   )
# 
# map

>>>>>>> ab384e6cd3f017dce9e7d273fac41bf48f398611